#---
#FUNCTIONS: ESTIMATE INTENSITY FUNCTION USING EMS ALGORITHM
#---

#quad5pt: computes Gaussian quadrature for 5 points
quad5pt <- function(J)
{
	JintervalWidth <- diff(c(0, J))
	nJ <- length(J)

	absc <- c(-0.906179845938664,
        -0.538469310105683,
		0.000000000000000,
        0.538469310105683,
        0.906179845938664)
	wgt <- c(0.236926885056189,
        0.478628670499366,
        0.568888888888889,
        0.478628670499366,
        0.236926885056189)

	nQuadpts <- length(absc)
	halfw <- 0.5 * JintervalWidth
	midpts <- 0.5 * (J + c(0, J[-nJ]))
	x <- rep(midpts, each = nQuadpts) + rep(halfw, each = nQuadpts) * rep(absc, nJ)
	return(list(nQuadpts = nQuadpts, lattice = x, wgt = wgt))
}

#quad10pt: computes Gaussian quadrature for 10 points
quad10pt <- function(J)
{
	JintervalWidth <- diff(c(0, J))
	nJ <- length(J)

	absc <- c(-0.973906528517172,
        -0.865063366688985,
        -0.679409568299024,
        -0.433395394129247,
        -0.148874338981631,
        0.148874338981631,
        0.433395394129247,
        0.679409568299024,
        0.865063366688985,
        0.973906528517172)
	wgt <- c(0.066671344308688,
        0.149451349150581,
        0.219086362515982,
        0.269266719309996,
        0.295524224714753,
        0.295524224714753,
        0.269266719309996,
        0.219086362515982,
        0.149451349150581,
        0.066671344308688)

	nQuadpts <- length(absc)
	halfw <- 0.5 * JintervalWidth
	midpts <- 0.5 * (J + c(0, J[-nJ]))
	x <- rep(midpts, each = nQuadpts) + rep(halfw, each = nQuadpts) * rep(absc, nJ)
	return(list(nQuadpts = nQuadpts, lattice = x, wgt = wgt))
}

#ems.constructor: computes smoothing matrix
#panel.data - panel data
#bandwidth - bandwidth
#quad.order - number of points
ems.constructor <- function(panel.data, bandwidth, quad.order = NULL, J = NULL)
{
	nInterval<- nrow(panel.data)
	if(is.null(J))
	{
		J <- sort(union(panel.data$L, panel.data$R))
	}

	if(any(J==0))
	{
		J <- setdiff(J,0)
	}
	
	if(quad.order == 10)
	{
		quadRule <- quad10pt(J[which(J!=Inf)])
	}
	
	else
	{
		quadRule <- quad5pt(J[which(J!=Inf)])
	}

	#construct the incident matrix
	censor.mat <- outer(tapply(panel.data$CT, panel.data$ID, unique), J, ">=")
	offsets <- apply(censor.mat, 2, sum)
	tmp.data <- panel.data[which(panel.data$M!=0),]
	incident.mat <- outer(tmp.data$R, J, ">=") - outer(tmp.data$L, J, ">=")
	  
	#construct the smoothing matrix at different x (for the final intensity estimate)
	tmpJ <- J[which(J!=Inf)]
	JintervalWidth <- diff(c(0,tmpJ))
	lattice <- quadRule$lattice
	tmpK <- outer(tmpJ, lattice, function(x, y) pnorm(q = (x - y)/bandwidth)) - 
		outer(c(0,tmpJ[-length(tmpJ)]), lattice, function(x, y) pnorm(q = (x - y)/bandwidth))
	K_x.mat <- diag(1/JintervalWidth) %*% tmpK
	if(any(panel.data$R == Inf))
	{
		K_x.mat<- rbind(K_x.mat, 0)
	}

	#construct the smoothing matrix used in the proposed EMS algorithm by Gaussian quadrature rule.
	wgt.mat <- kronecker(diag(0.5 * JintervalWidth), quadRule$wgt)
	K.mat <- K_x.mat[,(1:length(quadRule$lattice))] %*% kronecker(diag(.5*JintervalWidth), quadRule$wgt)

	if(any(panel.data$R == Inf))
	{
		tmpK <- 1 - t(outer(max(tmpJ), quadRule$lattice, function(x, y) pnorm(q = (x - y)/bandwidth)) %*% kronecker(diag(.5, nrow=length(tmpJ)), quadRule$wgt))
		K.mat <- cbind(K.mat, c(tmpK,1))
	}

	#returns smoothing matrix
	return(list(J = J, incident.mat = incident.mat, offsets = offsets, lattice = lattice, K_x.mat = K_x.mat, K.mat = K.mat))
}

#ems.nhpp: iteratively compute intensity function estimates for nonhomogeneous Poisson process
#panel.data - panel data
#bw - bandwidth
#quad.order - number of points
ems.nhpp <- function(panel.data, bw, tol = 1e-12, quad.order = NULL, partition = NULL)
{
	ems.input <- ems.constructor(panel.data, bandwidth=bw, quad.order=quad.order, J=partition)
	offsets <- ems.input$offsets
	tmp.data <- panel.data[which(panel.data$M != 0),]

	#EM iteration
	Lambda.old <- tmp.data$M %*% ems.input$incident.mat %*% diag(1/offsets)
	Lambda <- matrix(0, nrow=1, ncol= length(ems.input$J))
	niter <- 1
	repeat 
	{
		V.inv <- diag(1/as.vector(Lambda.old %*% t(ems.input$incident.mat)))
		tmpEM <- tmp.data$M %*% V.inv %*% ems.input$incident.mat %*% diag(as.vector(Lambda.old)) %*% diag(1/offsets)
		Lambda <- tmpEM %*% ems.input$K.mat

		if(sum((Lambda.old - Lambda)^2) > tol)
		{
			Lambda.old <- Lambda
			niter <- 1 + niter
		}
		
		else
		{
			lambda <- tmpEM %*% ems.input$K_x.mat
			break
		}
	}

	return(list(Lambda.est = as.vector(Lambda), lambda.est = as.vector(lambda), nIter = niter, lattice = ems.input$lattice))
}
